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Using an interpolant form for the gradient of a function of position, we write an integral version of 
the conservation equations for a fluid. In the appropriate limit, these become the usual conservation 
laws of mass, momentum and energy. We also discuss the special cases of the Navier-Stokes equa- 
tions for viscous flow and the Fourier law for thermal conduction in the presence of hydrodynamic 
fluctuations. By means of a discretization procedure, we show how these equations can give rise 
to the so-called "particle dynamics" of Smoothed Particle Hydrodynamics and Dissipative Particle 
Dynamics. 
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<U ; I. INTRODUCTION 

Q 

. The inherent difRculties of the equations for fluid hydrodynamics has given rise to a variety of schemes that 
numerically simulate the behavior of a fluid at hydrodynamic scales Of renewed interest, there have appeared 
the simulation schemes based on the idea of substituting the fluid by "particles" that move under the influence of 
forces that, in a coarse-grain limit, reduce to some approximate form of the hydrodynamic Navier-Stokes equations; 
these schemes are called Smoothed Particle Hydrodynamics (SPH) and Dissipative Particle Dynamics (DPD) 
Q. Although their numerical implementation is somewhat different, the latter including a random force, we shall 
\ argue that their origin is essentially the same. The use of these particle-like simulations has been reported of being 
P\| i: successful in different applications of fluid dynamics . An attractive feature of these simulations is that 

T-H ' one can use the enormous experience gained from the techniques of simulations of standars molecular dynamics; in 
C ■ particular, in dealing with rheological fluids. 

In this article, we present an integral representation of the hydrodynamic conservation laws based on the concept 
of the interpolant of a function |^,^ ; the interpolant is an integral representation of a function in terms of a weighting 
function. This function, in the appropriate limit, becomes a delta function and the interpolant yields an identity. 
C _ For calculational purposes one does not take such a limit and, therefore, one ends with an approximating form of the 
. corresponding function. The SPH simulation is based on a certain form of interpolant of the hydrodynamic functions, 
(~| ' such as density, velocity, etc. Here, we base our scheme not on the interpolants of the functions, but rather on the 
O interpolants of the gradient of the functions. This minor change proves to be very useful in writing integral equations 
O that reduce, exactly, to the conservation laws in the limit. Moreover, we shall see that the conservation of mass, 
^ , momentum and energy is exactly preserved at the level of the integral forms. At an approximation level, we shall 
show that the integral forms of the conserving currents are correct up to second order terms in the gradients of the 
corresponding fields. Since the integral expressions can be written down following the forms of the true laws, the 
■ phenomenological variables, such as viscosities and thermal conductivity, can be naturally included. In the same 
" " ' fashion, the extension to fluctuating hydrodynamics can be readily performed. 

The resulting integral conservation laws may then be used as an alternative to the exact differential laws. Further, 
they can be utilized as the starting point for approximate solutions. In this context, we show that by an appropriate 
discretization of the integrals one can render the equations to look as the equations, not for the hydrodynamic fields 
evaluated at space-time points, but rather, for quantities pertaining to particles. For instance the field velocity v(r, t) 
becomes the velocity Vi(i) of a particle at position Vi{t); a law of motion for the latter must then be supplied. 
Within the present theory one readily finds the equations of motion for the particles that are fully consistent with the 
hydrodynamic equations. We shall show how the SPH and DPD equations may be found. As we shall see, in general, 
one can have additional terms arising from the convective non-linear terms of the hydrodynamic time derivative. 

We organize the article as follows. In the next section we introduce the interpolant of the gradient of a function. 
With such an object, we write down the conservation laws, both in general and for the cases of Navier-Stokes and 
Fourier law. The extension to fiuctuating hydrodynamics is also shown. In section HI we present a discretization 
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procedure that yields a particle-like simulation algorithm. We discuss the relationship with the SPH and DPD theories. 
We conclude and make additional remarks in section IV. 



II. HYDRODYNAMICS IN INTEGRAL FORM 



A. An interpolant of the gradient 

The idea of formulating an integral version of the laws of hydrodynamics is based on a limiting representation of the 
gradient of a function of position r. This representation we call it an interpolant, following Lucy j^] and Monaghan 
1^ SPH formulation. As we mentioned in the Introduction, we use the interpolant of the gradient of a function rather 
than that of the function itself. 

First, we show that the following identity is correct, 

\7A{r) = -lim,„_o J dr' A(r') I^-lJ|vF(|r - r'|; ro), (1) 

where the tensorial character of A is left unspecified, and M^(|r — r'|;ro) is a distribution or weighting function sharply 
peaked with width rp. We demand that all its moments exist although the function itself may not be integrable. We 
assume the following: 

" dr\rrW{\r\;ro)^r^~^Mn (2) 

with n > 1. For n ~ 1 we require Mi = 3 but for n > 1 we leave A4n unspecified. These requirements are easily 
satisfied by noticing that the integrand of (|l|) can be related to the gradient of a distribution that tends to a delta 
function, namely 

^iy(|r|;ro) = -V/(|r|;ro) (3) 

where /(r; tq) is such, that in the limit tq ^ becomes 

\imr„^of{r;ro) ^ ^S{r). (4) 

Clearly, any distribution that tends to S{x) can be used. It is important to stress that VF is a function of the magnitude 
r| and not of the corresponding vector. This dependence is important for the use of approximations and, as we shall 
see later on, for the setting up of the conservation laws. 

For the validity of (|l]) we further require that A{r) is analytic everywhere inside the domain of the integral. This is 
not a stringent requirement since we are interested in hydrodynamic fields. Calling I(A(r)) the integral in equation 
(|l|), we first make a change of variable, r' ^ r + r', and then, we perform a Taylor expansion of A{r + r') around r. 
We obtain 

/(A(r)) = 5,A(r) J dr' "^^W {\r' \; ro) + 

ia„a^9,A(r) J dr'^ii^|^W^(|r'|;ro) + ..., (5) 

where = d/dr^ with the cartesian components of the vector r. Due to the spherical symmetry of W all the odd 
powers of r' (even powers in the derivatives) vanish identically. With the use of Eq. (0) the first integral in the rhs of (|^) 
gives the Kronecker delta Sa/s, independently of tq, while the second integral yields rQ{Saf:!Sjri + Sa'ySf3r] + Sar)Sf3j)W3/15. 
Clearly the higher order terms are proportional to Tq" times odd (2n + l)-derivatives of A. We write, generically, 

I{A{r)) = Vyl(r) + O (r^V^A) . (6) 

In the limit rg — > 0, the interpolant is the gradient of A. For approximation purposes, we note from (^) that the 
correction to the gradient is third order in the derivatives. 
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B. Conservation laws 



The study of the hydrodynamics of a fluid is based on the conservation laws for mass, momentum an energy. 
In the following we assume that all the fields are evaluated at a spatial point r and that all are time dependent as 
well. The conservation laws are 



dp 
dt 
di 

dt 
de 

'dt 



-V • (pv) , 

-V • (^jv + PI - n 

-V • (ev + (PI - n) • v + J 



(7) 
(8) 
(9) 



where p is the mass density; v is the velocity of the fluid; j is the momentum density, j = pv; P is the hydrostatic 
pressure; 11 is the viscous stress tensor; e = + m is the total energy density with u the internal energy density 

of the fluid; and J is the heat current. Integration over the whole volume of the fluid shows that the total mass, 
momentum and energy of the fluid are conserved. 

With the use of the interpolant given by Eq. ([l]), we can write down analogous conservation expressions (with 
implicit time-dependence) : 



dp{v) 
dt 

dt 



de{r) 

~dr 



dr'll/(|r-r'|;ro)- 

dr'iy(|r-r'|;ro)- 
[P(r; r') + P(r';r) 

dr'VK(|r-r'|;ro)- 

+ [P(r;r')«a(r)+P 
+ [J,(r; r') + Ja(r';r)]}. 



r - r')c 



[/9(r')«a(r) +p(r)«„(r')] 



r')c. 



fp{[jXr')«a(r)+jMr)z;„(r')] 



|r - r' 

Sail - [n„0(r;r') +n„/3(r';r)]} 
r — r') 



(10) 



(11) 



^^{[e(r')«.(r)+e(r)«„(r')] 
r';r)t;„(r')] - [Iiap{v-v')vfi{v) +I{ap{v' ■v)vp{v')] 



(12) 



The above expressions, by the rule of the interpolant, are equal (in the limit tq 0) to (minus) the divergence of 
the terms in square brackets evaluated at r' = r; in particular, the "kernels" P(r;r'), ll(r;r') and J(r;r') must be 
chosen such that the interpolant equal the divergence of the actual pressure, viscous stress tensor and heat current 
when r' = r. But before we discuss how to choose these kernels we point out that all the terms inside the square 
brackets are symmetric with respect to the interchange of the variables r and r'. Therefore, integration with respect 
to r makes the rhs of all the equations vanish, thus yielding the conservation of the extensive variables, independently 
of whether the limit ro — > is taken or not. 



C. Constitutive relations 



The conservation laws must be provided with constitutive relations in order to have a closed set of equations. For 
discussion purposes we shall choose the mass density, velocity and temperature r(r, t) as the independent fields. 
Therefore, for the pressure and the internal energy we need to know the equations of state of the fluid in terms of p 
and P; in particular, we assume we know the functional dependence P{p,T). Thus, the kernel for the pressure may 
be chosen as 

P(r;r')=P(p(r),P(r')). (13) 
Clearly, the kernel is not symmetric in its variables and hence we need the symmetrization in the conservation 



equations (10) - (|12|) above. Substitution of this form into, say, Eq. (|l l|) , gives to lowest order, the term 



dr' [P(rir') + P(r'ir)] jj--l2n/(|r - r'|;r„) = (14) 
'^^'vp-flSlvT, ,15, 



dp J \dT 
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where the rhs is evaluated at r' = r. The rhs is — VP(r). 

Regarding the viscous stress tensor 11, we use the usual one that gives rise to the Navier-Stokes equations, linear in 
the velocity gradients. The corresponding kernel may be written as 



n„0(r; r') - -r;(r) J dr"Wi\r' - r"|; m) 



(r' - r")a 



r' - r" 



C(r) / dr"M/(|r'~r"|;ro) 



r yr r )i, II 

'5"/3-F3 — zirrMr ) 



(16) 



where the viscosities 77(r) and C(r) are, either, given functions of r, or they depend on r through a further dependence 
on density and temperature. Substitution of this equation into Eq. (pT|) , yields to lowest order the familiar viscosity 
terms of the Navier-Stokes equations: 



dv' [n„/3(r;r')+n„/3(r';r)] 



l^W{\v-v'\-r,) 



ri{r) [ daVf}{r) + dfjVair) - -Sapd^v^{r) ) + C{r)da[3di,v^{r) 



(17) 



For the form of heat current J we consider Fourier law in which the current is linear in the temperature gradient. 
The kernel can be written as, 



J(r;r')=«(r) / dr"W^(|r' - r"|; ro) j^J— |JJ|T(r"). 



(18) 



Again, substitution into Eq. (h2|) yields, to lowest order. 



dr' [J,(r; r') + J,(r'; r)] • ^^^^W{\r - v'\;ro) = V • [^(r)Vr(r)] 



(19) 



where the r-dependence of k may be given through its dependence on temperature and density. 

Thus, we have shown a consistent way of presenting an integral form of the equations of hydrodynamics, such that, 
in the appropriate limit yield the true equations; we note that all the phenomenological coefficients are readily and 
unambiguously identified. This is an important point since, as we shall show in the next section, the choice of the 
functional forms of the stress tensor and the heat current are, by no means, unique; that is, we can prescribe different 
functional forms that in the limit also give similar expressions to the usual hydrodynamic laws. We shall defer further 
discussion of this point to section III. 



D. Hydrodynamic fluctuations 

We now turn our attention to the formulation of the study of hydrodynamic fluctuations. Following Landau and 
Lifshitz [ p^p^ we limit ourselves to small fluctuations around a given flow, solution to Eqs. (|^) - (^), so that a 
linearization in the fluctuations is possible. In keeping with our assumption that the indepedent variables are the 
mass density, velocity and temperature of the fluid, we define the fluctuations as linear deviations from the flow, that 
is 

p{r,t) ^ po{r,t)+dp{r,t) (20) 

and analogous expressions for v(r, i) and T{r,t). The functions po{r,t), Vo(r,t) and To(r, t) constitute a given flow, 
solution to the full non-linear integral equations ( p^ - ( p^ ) , with the expressions (|l5|) and (^8^ . Now, using expressions 
such as (po| ) we can linearize the integral equations in the fluctuations; for instance, the continuity equation for the 
fluctuations, cf. Eq. (p^) , becomes 



J dr'W^(|r-r'|;ro)|^ 



at J " " |r-r'| 

X [po(r')'5v(r) + (5p(r')vo(r) + po(r)<5v(r') + <5p(r)vo(r')] (21) 
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and similar linear equations for the (partial) time derivatives of momentum density and energy, 5j and 5e. Next, 
one identifies the source of the fluctuations as arising from spontaneous fluctuations of the stress tensor and the heat 
current. This is implemented in the usual way, by adding to the linearized equations for the fluctuations of the 
momentum and energy densities, terms proportional to the divergence of a random stress tensor and to the divergence 
of a random heat current, respectively. 

To be precise, we add to the equation for the fluctuation of the momentum density, 

/ dr'W^(|r-r'|;ro)^^— 47-fl«(r,r',t) (22) 
J |r-r'| 

where the tensor n-'^(r, r', t) is a gaussian random stochastic function, symmetric under interchange of r and r', with 
zero mean and with its second moment obeying the usual fluctuation-dissipation relations, 

(nf^(ri,r2,0n^,(r3,r4,O) = 
2fc{[To(ri)?7o(r2) + To(r2)??o(ri)] {Sa-fSpi, + 5ai,6p^) 

+fc[To(ri)(Co(r2) - ^^(ra)) + To(r2)(Co(ri) - ^?7o(ri))]<5„/3^.} 

x6{t - t') {Sin - r3)(5(r2 - r4) + S{ri - ri)d{r2 - rg)} . (23) 

In the same fashion, we add to the equation for the fluctuation of the energy density, 

/ dr'W{\r-r'\;ro)p^-3''ir,r',t) (24) 
J |r-r'| 

where J^(r, r', t) is a gaussian random stochastic function, symmetric under interchange of r and r', with zero mean 
and with second moment 

(Ji^(ri,r2,0J#(r3,r4,t')) = 
2k (r2(ri)Ko(r2) + T2(r2)Ko(ri)) Sa^pSit - t') 
X [Sin - r3)(5(r2 - r4) + S{ri - r^)S{v2 - rg)] . (25) 

We recall that the quantities Sp{r,t), 5j{r,t), etc. depend from those of the underlying flow, po{r,t), jo(r,t), etc. but 
not the other way around. That is, one first solve for the latter and then one finds the fluctuations. It is understood 
that the flow is stable; that is, one should always have Sp/po <C 1, etc. In the next section we shall also comment how 
one can include the fluctuations within a particle-like simulation. 



III. A DISCRETIZED INTEGRAL HYDRODYNAMICS 



The integral formulation presented in the previous section is simply an approximate representation of the usual 
hydrodynamic laws. Their usefulness resides on whether their solution may be easier to find than those of the actual 
equations or on their amenability for approximations. In this regard, we recall the approximate SPH and DPD schemes 
where a particle-like simulation, similar to a molecular dynamics simulation, represents the flow of a continuous fluid. 
In this section we present a particular discretization of the integral equations of section II that may be used as the 
basis for a simulation in terms of "fluid particles." 



A. A particle-like scheme 

The basic idea is first to divide space in cells of finite size AV^ and, then, to define the field variables for each cell. 
We call Vi the position vector of the i-th cell, and the following list summarizes the variables for such a cell: 

p(r,t)AV mi{t) mass 

j{r,t)AV Pi{t) momentum , , 

e{r,t)AV <ii{t) energy ^ ' 

v(r,t)Ay Vj(i) velocity 
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Now, the kernels of pressure and viscous stress become "potentials" of force between the i-th and j-th cells while 
the heat density current is now a "current" of energy between such cells: 



P{r;r';t)AV 
[I"^(r;r';t)AV 
J^{r;r';t)AV 



Vij (t) pressure potential 
U°'>^{r;r';t)AV KfW stress potential 



(27) 



heat current 



The integrals are then discretized by summing over cells directly and not over labels that localize the cell in a 
Cartesian grid, that is. 



N 

J dr^^AV^, 



(28) 



where we have assumed there are N cells in the total volume. This discretization implies a careful choice of the 
discretized version of the weighting function W{\r — r'|). That is, we cannot simply change r by and r' by Vj in 
the functional form of W since the equations for the moments, Eq. (^, would not be correct. This is due to the fact 
that those results make use of the spherical symmetry of W. Instead, we propose the following discretization that 
give rise to the correct moments: 



W{\r - r'\)AV >V(ry) = ATrrfJAVy/^ W{nj) 



(29) 



where we have defined rij — jr; — rj\. This form also takes into account that W is always part of an integrand. As a 
particular example, using as a representation of a delta function. 



Sir) 



lim,„^o-e""/''°, 
ro 



yields for W{r), 



and, correspondingly for W{i 



^3 >^ 



W{r) 



W(r,,) 



47rr4 



{AV) 



1/3 



1 



which shows that, in discretized form, all the moments but the zeroth are well defined. 
With the above reformulation, and defining 



(30) 



(31) 



(32) 



(33) 



the conservation equations look as follows. Conservation of mass: 



(34) 



Conservation of momentum: 



Conservation of energy: 



dpf 



dej 
'dt 



j 

- {^tj ■ V.» + • Vj) + {J,j + Jji)} . 



(35) 



(36) 
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By construction, the total mass, total momentum and total energy are conserved. This can be seen by summing 
the above expressions over i = 1, 2, . . . , TV. 

So far, the equations are quite general and one needs constitutive relations for the kernels of pressure, viscous stress 
and heat current. For instance, the viscous stress tensor linear in the velocity gradients may be found by discretizing 



Eq.®, 



k 

k 



a , a ^ ca0 v 



C.5"^VW(r,fe)e,-,«,^ (37) 



while the heat current may be found from Eq.(p^) 

J^^n,^W{r,i,Ye,uTk. (38) 
fe 

One should keep in mind that in order to close the equations one still needs the equations of state for the pressure 
Vij = V(mi, Tj) and the internal energy per particle = u{rai, Tj). 

Up to here it is simply a discretization of the equations. The interesting addition now [§j3[^ is to assume that 
the positions of the cells become the positions of particles that are allowed to move. This is certainly a bold 
assumption, since making the fluid particles move should be done in a Lagrangian formulation of the fluid dynamics 
rather than in an Eulerian one. The present are a different type of particles, however, since as we can see from the 
above equations, they not only change their momenta but they also have variable mass and carry with them their 
internal energy in addition to their kinetic one. Nevertheless, one may justify it by arguing that one is actually looking 
at every instant of time to a state of the fluid not on a grid but rather on a "fluidized grid" ; the motion law of 
being used as an updating of the state of the fluid on such a grid. In any case one can asses the validity of such an 
assumption a posteriori; as we mentioned in the Introduction, successful simulations of actual flows with schemes like 
the present one have been reported, see Rcfs. j5|J^,[7|,p|,p|jlC|] . Although the choice of the law of motion is arbitrary, it 
seems "natural" to consider the rate of change of the position of the cell as the velocity of the fluid at that point: 

—^=-^{t). (39) 

We point out that this choice is not unique; see Ref. |^ for other forms used in SPH simulations. 

In principle the above scheme is complete and closed. However, for an actual implementation of a simulation 
based on it there are further questions to resolve, such as the boundary conditions in terms of the particles and 
the discretization of time. Since there are already in the literature a host of procedures ||7|,p^ both for dealing with 
boundary conditions between particles and solid frontiers and for the time discretization, we shall only discuss the 
latter because of its relevance in the inclusion of hydrodynamic fluctuations. 

A simple algorithm to simulate the dynamics consists of a two-step propagation in time Q]. First, there is a 
"collision" step in which one finds the values of rui, Pi and at time t + At from from knowledge of all the variables 
at time t, using Eqs. (|3) - (|3|), with 



dA,[t) A,{t + M) - A,{t) 



(40) 



dt At 

It is then followed by a "propagation" step in which the positions {t + At) are computed using 

r, {t + At) « r, (t) + At (t + At) . (41) 

This combination is more accurate than if both steps were done with Eq. (^0|). This algorithm, however, is also 
useful to include the fluctuations as part of the evolution and not as a posteriori calculation, being thus helpful in 
determining the stability of the flow. This is an important point since the purpose of the simulations is to solve the 
equations by an actual propagation in time of the flow. (An analytical solution need not be done in this way; for 
instance, if the equations are linearized one may solve them using an integral transform technique.) Thus, one can 
include in the equations for the momentum and the energy, Eqs. ( ^5| ) and (|36|), discretized versions of the random 
viscous tensor tt,^, and random heat current J^^ , both symmetric in ij. Since their second moments must obey 
discretized versions of the fluctuation-dissipation expressions (psh and (p5|) these tensors may be added as 
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n^^-^it) = {2m{tH{t)+2kT,{tH{t)f'AXrfit) 

+ {mt)Q{t) + kT,{t)Ut)f^ ^Yf{t) (42) 

(0 = im{t)^,{t) + kT,{t)K.,{t)f/^ ^Z^,{t) (43) 

where AX"'^^ (t), AY^"^{t) and AZ°j{t), symmetric in ij, represent independent random increments (Wienner processes 
) with zero mean and correlations 

l^AYf{t)AYZ{t')) = 5''^ 5'"' {SaSjm + ^.™<5,/) 5u' 

(^AZf^{t)AZ^.^{t)) = 5'-^ {6a5j„. + 5^.^5,l) Sw (44) 

From a practical point of view, any good commercial pseudo-random number generator suffices for these increments. 
The important aspect we want to stress is that, in contrast to the continuous version (|2^) and (p5|), here the tem- 
perature, viscosities and thermal conductivity that appear in Eqs.(^2|) and ( ^3|) are evaluated at the current values 
of the full fluctuating quantities and not at the values of the variables at the underlying flow. Since it is assumed 
that the fluctuation are always small and do not make the flow unstable, this is a minor approximation. Moreover, 
as mentioned above, if the flow do become unstable by adding the random viscous tensor and heat current that 
may imply, barring numerical inaccuracies, that either the flow is indeed unstable or that the method itself does not 
faithfully describe the flow. 

It should be clearly understood that the particle-like representation of a continuous fluid flow depends on two 
different approximation; first, one approximates the true differential laws by integral expressions with a finite width 
To of the weighting function; and second, the integrals are discretized. These approximations pose constraints on 
the length scales of the fluid. On the one hand, the density of point-particles must be such that the mean particle 
separation {AV)^/'^ is smaller than tq in order to have a good approximation of the integrals; and on the other hand, 
a typical hydrodynamical length, call it A, must be larger than tq itself in order to have a good representation of the 
gradients in terms of the integrals (i.e. an independence of the parameter tq, see Eq.(|6|).) That is, one should always 
have, 

(Al/)i/3 < ro < A. (45) 

This way, the limit tq ^ implies not only the equality of the integral and differential forms of the conservation laws, 
but also the continuum limit itself. 



B. SPH and DPD as special cases 

The purpose of this section is not to make a revision nor a comparison of SPH and DPD schemes with the present 
one, but rather, to show that they may be viewed as special cases of a more general scheme that reduces to the 
macroscopic conservation laws of fluids. 



As we have seen the discretized conservation equations (|3j)-(36) are very general. One still needs to provide 
constitutive relations for the pressure, viscous tensor and heat current and, as long as ij-symmetrized forms are 
provided, the conservation laws are guaranteed. The particular expressions given in the previous section, such as ( |3^ ) 
and (^Sf), are just examples. But before we present other forms used, such as those of DPD and SPH, we want to 
mention some aspects of the time derivatives used. 

In the schemes used in SPH and DPD, the time derivatives of the properties of the particles have been interpreted 
as already including the convective contribution. In the discretized version this is equivalent to identify 

^^>V(ry)ey -v,™, (46) 

j 

5]W(r.,)ey--v,p« (47) 



dmi 


dmi 


IT 




dpl 


dp? 


dt 


dt 


dei 


dfi 


'dt 


^ dt ' 





3 



W(r,j)e,j • v,ej (48) 
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This is a subtle point: One the one hand, one could argue that the derivative is following the motion of the fluid 
particle, like in a Lagrangian formulation; however, the right-hand-side of the corresponding conservation laws (0)- 
( ^ ) should be accordingly transformed. Since the latter is not done in SPH and DPD, one may still say that those 
formulations correspond to not too large Reynolds numbers where one can neglect the convective contributions. It 
may be interesting to include those terms explicitely in a simulation. 

With the above identification of the time derivatives, we can now see a closer resemblance to the equations of SPH 
and DPD. For purpose of exemplifying the relationship, we shall only discuss the equation for the momentum. Using 
the mass conservation equation and the fact that = rriivf, the equation for the momentum can be written as, 

= E ^('^H {^"' (^^. + ^0 - (4' + (49) 
j 

In Refs. SPH is formulated with forms for the pressure such as 

^,+7',,-^^ + ^^ (50) 

^ J J '■ 

with a and h constants and with a given equation of state for Pi in terms of rrii . jl5| 
The viscous stress tensor of SPH and DPD may be generally written 

+ nff = A(r,,)e^(vf - v^) + i?(r,,)e§« - vf ) (51) 

with appropriate choices of A and B [^PJ^,P|,|To| . It is interesting to note that this form can give rise, in the continuum 
limit ro — > 0, to terms proportional to V^v and V(V-v); however, one cannot independently identify the corresponding 
viscosity coefficients. This is to be contrasted with the expression of the tensor given by equation (|3^) where there is 
an independent identification of the viscosities. 

In the DPD simulations there is an additional ingredient. Namely, that the pressure term is taken to be stochastic. 
Within the present scheme this may be interpreted as including hydrodynamic fluctuations with white noise and with 
a particular temperature and viscosity as given by equation (p3). In this regard we differ from the interpretation of 
DPD equations given by Espahol et al. ||^ and Marsh et al. |lO[|. In that interpretation, the equations of DPD are 
taken as the "microscopic" dynamics of the particles of a fluid, from which the macroscopic laws are to be extracted, 
much in the spirit of Langevin and Boltzmann equations. Within that interpretation they have argued that the 
random part should be modifled in order to account for the correct fluctuation-dissipation relation of Langevin-like 
equations. According to the present theory, since the equations of motion of the particles are only an approximation 
to the macroscopic equations of the fluid flow, there is no need to modify the DPD equations. Therefore, the random 
contributions of DPD may be seen to already refer to hydrodynamic fluctuations. Moreover, if one wishes to find 
the corresponding Fokker-Planck equations to the discretized hydrodynamic equations ( p^ - (^6|) one can follow the 
theory of van Saarloos et al. ||l^ of non-linear hydrodynamic fluctuations. 



IV. FINAL REMARKS 



In this article we have presented an integral form of the conservation laws of a macroscopic classical fluid in terms 
of an interpolant for the gradient of a given function of space. This form is amenable for a discretization of space 
and may be interpreted in terms of the dynamics of "fluid particles" . To complete this discretized dynamics one must 
provide the law of motion of the position of the particles; one may prescribe that the field velocity equals the rate of 
change of the position of the particle. Within this scheme one can easily find out the corresponding Navier-Stokes 
equations of viscous fiow and the Fourier law of heat conduction. Moreover, hydrodynamic fiuctuations can also be 
readily taken into account. 

We have argued that numerical simulations currently used, known by Smoothed Particle Hydrodynamics (SPH) 
and Dissipative Particle Dynamics (DPD) 0|, and which are based on a particle- like simulation of a continuum fiuid, 
may be seen as special cases of the general formalism here presented. In that way, one the one hand, one guarantees 
that the simulations have a correct continuum limit, and on other, there is a clear route of how to represent, in the 
particle dynamics, known effects of macroscopic fluids; for instance, with the present theory one can see how to include 
thermal effects, absent in DPD simulations. 

Finally, we want to stress the potential uses of this type of schemes. It does seem that a complication in any 
simulation of fluids is the discretization of space with its concomitant difficulties of boundary conditions; this is the 
more important if one is interested in rheological fluids, such as suspensions. That is, in order to simulate a simple 
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flow a discretized differential scheme (e.g. finite difi^erences) may appear to be better than a discretized integral 
version; this is because the latter makes use of the weighting function which in turn must resemble a delta function, 
and therefore, it appears that one needs more "particles" than points in a grid Q, cf. Eq. (^). However, having paid 
this price, there is a host of "tricks" and techniques, borrowed from standard molecular dynamics, that can be used 
to simulate moving boundaries and solid objects, e.g. Lee and Edwards shear boundaries |T^ ], or "freezing" certain 
number of particles to simulate a rigid body etc. Moreover, the inclusion of hydrodynamic fluctuations also seem 
to be much easier within a particle-like scheme than within a field-like. 
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